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Abstract 

Previous numerical studies have revealed the existence of embedded solitons (ESs) in a 
class of multi-wave systems with quadratic nonlinearity, families of which seem to emerge 
from a critical point in the parameter space, where the zero solution has a fourfold zero 
eigenvalue. In this paper, the existence of such solutions is studied in a three-wave model. 
An appropriate rescaling casts the system in a normal form, which is universal for models 
supporting ESs through quadratic nonlinearities. The normal-form system contains a single 
irreducible parameter e, and is tantamount to the basic model of type-I second-harmonic 
generation. An analytical approximation of WKB type yields an asymptotic formula for the 
distribution of discrete values of e at which the ESs exist. Comparison with numerical results 
shows that the asymptotic formula yields an exact value of the scaling index, —6/5, and a 
fairly good approximation for the numerical factor in front of the scaling term. 
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Soliton solutions in nonlinear models are parameterized by their intrin- 
sic frequency. An elementary principle determines intervals of values of the 
frequency where ordinary solitons may be found: these must lie in gaps 
(semi-infinite or finite ones) in which the linearized version of the corre- 
sponding model has no solutions for radiation waves. A usual argument 
supporting this principle is that if one tries to construct a soliton with the 
frequency falling within a band occupied by the radiation-wave solutions, 
the soliton would start to decay by emitting linear waves whose frequency 
is exactly equal to the soliton's intrinsic frequency. Nevertheless, it has 
been demonstrated that exceptions are possible in practically important 
models (notably, in a system of three waves coupled by quadratic nonlin- 
earity, with an extra linear coupling between two fundamental-frequency 
waves induced by a Bragg grating, in an underlying model of the spatial 
evolution of optical fields in a planar waveguide). These exceptional soli- 
tons, which are embedded in the radiation spectrum, exist as a series of 
isolated solutions at discrete values of the frequency inside the radiation 
bands. An issue of considerable interest is the evolution of these solutions 
following variation of governing parameters of the model: in that sense, the 
embedded solitons form continuous families. Previously, it was found in an 
empirical form (from direct numerical results) that a presumably infinite 
series of embedded soliton solutions emanates from a single point in the 
parameter plane of the above-mentioned three-wave system. The nature of 
this phenomenon was not understood. 

In this work, we resolve these issues, by developing an asymptotic anal- 
ysis and verifying the predictions against numerical results. As a result, 
we derive a universal asymptotic approximation (normal form) for models 
with quadratic nonlinearity that support embedded solitons. The normal 
form amounts to a system of two second-order equations, which is well 
known as a basic model for ordinary solitons in second-harmonic-generating 
systems. There was a rather common belief that all soliton solutions had 
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been found in this much-studied system. Nevertheless, in this work we are 
able to predict the existence of an infinite series of previously unknown 
embedded solitons in the model. This is done in an analytical form, by 
realizing that the embedded solitons feature a broad inner zone, the so- 
lution in which must be matched to exponentially decaying ones in outer 
zones. Further, a method (Bohr-Sommerfeld quantization rule) borrowed 
from quantum mechanics is applied to the solution in the inner zone. The 
most essential part of the eventual analytical result is an asymptotic dis- 
tribution law for values of the single control parameter of the normal-form 
system at which the embedded solitons exist. Comparison with numerical 
results clearly shows that the analysis has produced an exact value of the 
respective scaling index. 

I. INTRODUCTION 

This article deals with spatial solitary waves (or solitons in short, without implying 
integrability of the underlying model) in a model of a planar waveguide with a quadratic 
(x 2 ) nonlinearity, where two fundamental-frequency (FF) waves, iiiand 11,2, are linearly 
coupled by the Bragg reflection from a grating of ridges or grooves on the surface of the 
waveguide. The waves are coupled into the waveguide so that their amplitudes evolve 
along the direction z, which is aligned with the Bragg grating, the Poynting vectors of 
the two FF waves having angles ±a with the z axis (where the angle a is small). The 
;\;( 2 )-induced nonlinear coupling between the FF waves generates a second-harmonic 
(SH) wave, u 3 . The corresponding dimensionless model, derived in Ref. amounts 
to a system of normalized equations, 



2i(u 3 ) z - qu 3 + D(u 3 ) xx + u x u 2 = 0, 

where x is the transverse coordinate, the opposite signs in front of the walk-off terms 
(^1,2)3. correspond to the above-mentioned opposite angles ±a (a itself does now appear 



i(ui)z + i(ui) x + u 2 + u 3 u* 2 = 0, 
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in the equations as it can be removed by a rescaling), the asterisk stands for complex 
conjugation, and real coefficients D and q account for the intrinsic diffraction and wave- 
number mismatch of the SH wave (the intrinsic diffraction of the FF waves is neglected, 
as the artificial diffraction, induced through the interplay of the Bragg-grating-induced 
linear coupling and opposite walk-off terms, is stronger by several orders of magnitude). 
Note that is invariant with respect to both a shift of phase by an arbitrary real 
constant 0, 

S : (ui,u 2 ,u 3 ) i-> (exp(i0)tii,exp(i0)w 2 ,exp(220)w 3 ), (2) 

and to an arbitrary transverse shift. 

Solitary waves in Eqs. (0) have been investigated in several works Q,Q|. The 
present work was motivated by results in the two latter references, which had produced 
embedded solitons (ESs) in the system. ESs are distinguished by the feature that their 
internal frequency (spatial frequency or propagation constant, in the present context) 
is in resonance with radiation waves, \^ . In general, one expects to find quasi-solitons 
with non-vanishing oscillating tails in such a case. However, at certain discrete values 
of the frequency, the amplitude of these tails can exactly vanish, thus giving rise to a 
truly localized solution with a frequency that is embedded in the continuous spectrum. 

One may seek for stationary ES solutions of (JIJ in the form 

u h2 (x,z) = exp(ikz)Ui i2 {€), u 3 (x, z) = exp(2ikz)U 3 (£) (3) 

where k denotes the wave number, and £ = x — cz is the tilted ("walking") coordinate, 
with c measuring the walk-off of the spatial soliton's axis relative to the propagation 
direction z. Substitution of the ansatz (JHJ) into (JIJ leads to a system of ordinary 
differential equations for the complex amplitudes L/j, 

-w x + z(i - C )u[ + u 2 + u 3 u; = 0, 

-kU 2 - z(l + c)U' 2 + U x + U 3 U{ = 0, (4) 
-(4Jfc + q)U 3 + DU'l - 2icU' s + U X U 2 = 0, 

with the prime standing for d/d^. 

The analysis in this paper will concentrate on standing solitons with c = 0, in which 
case we can set D — 1 with no loss of generality. Furthermore, Eqs. (j3J) with c = 
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FIG. 1: Enter Figure ^ here. 

allow for symmetry reductions to four- dimensional invariant subspaces, the phase-shift- 
invariance J2J) accounting for the existence of a one-parameter family of such invariant 
spaces. Localized solutions can be studied in any member of this family. A convenient 
reduction is achieved by setting U\ = —U2 = V = V\ + %V% (Vi and V2 are real), 
Us — C/3 = W (i.e., W is real), which reduces Eqs. (jlj) to a system 

-W + iV - V W - V* = 0, 

(5) 

-{Ak + q)W + W" - \V\ 2 = 0. 

A linear analysis reveals the region in parameter space where ESs may be found. 
Indeed, in order for ESs to exist, the zero solution of Eqs. (jJJJ), V = W = 0, has to 
be a saddle- center, i.e., the spectrum of the linearized system has to contain a pair 
of real eigenvalues, such that solutions exponentially localized in one direction exist, 
and a pair of purely imaginary eigenvalues, that account, in the general case, for non 
vanishing tails in another direction (if all the eigenvalues are real, the system gives 
rise to ordinary gap solitons instead of ESs, 0). For Eqs. ©, it is straightforward 
to find that the zero solution is a saddle-center in the region of the parameter space 
determined by the inequalities 

k 2 < 1, Ak + q < 0, (6) 

see also Fig. [T] 

An important property of Eqs. © is reversibility. This means its invariance under 
the transformation (involution) 

R X ;{V,W,W')»(V*,W,-W), (7) 

An ES of Eqs. ([TJ must be a solution of Eqs. (JHJ homoclinic to the saddle-center 0. For 
such a solution to exist, the one-dimensional stable and unstable manifolds of the equi- 
librium have to agree along the solution, and in general this situation is of codimension 
three in M 4 . However, in an irreversible system obeying the invariance (J7J), the stable 
and unstable manifolds of are i?i-images of one another. Thus, a homoclinic orbit to 

6 



Accumulation of embedded solitons 



exists if the unstable manifold of intersects the two-dimensional fixed subspace of 
the full four- dimensional phase space of Eqs. (JSJ), Fix(i?i) = {V 2 = W = 0}. This is a 
codimension-one situation and consequently, ESs should exist for parameter values on 
curves in the two-dimensional parameter plane of the system (|5|). Homoclinic orbits 
that intersect the fixed space of the involution (J7J) are themselves invariant under its 
action, therefore they are called symmetric. 

Remark 1 In addition, system is reversible with respect to a second involution, 

R 2 :(V,W,W')^(-V*,W,-W), £^-£. (8) 

According to the behavior of the first FF wave, V\, we will call ESs that are symmetric 
with respect to R\ even and those symmetric with respect to R 2 odd. (Note that the SH 
(W ) component of both R\- and R2-symmetric homoclinic solutions is automatically 
even.) 

Remark 2 This paper deals with fundamental ESs of (QJ), i.e., one-pulse (single- 
humped) homoclinic solutions of Eqs. ^ or It is well known that N -pulse homo- 
clinic orbits can emerge in bifurcations of homoclinic orbits to saddle-center equilibria 
. In particular, it can be shown that families of_2-_pulse homoclinic orbits (bound-state 
solitons) must accumulate on a family of ESs Jfl, \2j- 

It is straightforward to numerically construct both even and odd ESs, using a variant 
of a numerical shooting method proposed in jjj, which searches for intersections of the 
unstable manifold of with the corresponding fixed spaces Fix(i?i ) 2). Afterwards, 
branches of these solutions can been continued using the software "package AUTO M . 

nn 

The corresponding results have already been reported in Refs. |3|, ]^ ■ 

Figure El gives a concise overview of the computed branches of ESs for Eqs. (JHJ). In 
this paper, we are interested in the situation close to the line 4k + q = 0. If this curve 
is crossed (with \k\ < 1), the equilibrium changes its type from a saddle-center into 
a saddle. It has been shown in Ref. that there is a possibility for curves of ES 
solutions to approach this curve transversally. In that case, if the curve is crossed, the 
soliton is not destroyed but rather changes its type from embedded to an ordinary gap 
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FIG. 2: Enter Figure El here. 
FIG. 3: Enter Figure |3] here. 

soliton. As can be seen in FigureEl this scenario does not occur for system (JSJ). Instead, 
one finds that all the curves of the ES solutions approach the point k = l,q = —4, 
becoming asymptotically tangent to the line 4k + q = 0. For these parameter values 
is a degenerate co dimension- two point, whose linearization produces a fourfold zero 
eigenvalue. Our aim in this paper is to understand how families of ESs accumulate at 
(emanate from) this point. 

In order to understand this accumulation we will perform an asymptotic analysis 
of ES solutions. The analysis, however, will not be performed directly for Eqs. (j3J). 
Instead, we perform a rescaling of the original problem and derive a suitable normal 
form in Section |H] Afterwards, in Section UTTl an analytical WKB approximation yields 
asymptotic estimates for parameter values at which ES solutions are predicted to exist. 
The theory is verified by numerical results. 

II. THE TWO- WAVE SYSTEM AS A NORMAL FORM 

The analysis of Eqs. © is facilitated by a rescaling of the equations, suggested 
by the numerically computed solutions. First, we display examples of such solutions. 
Fig. El shows the FF and SH components of ES solutions, which have been computed 
along the line k = 0. It can be seen that the FF component features internal oscilla- 
tions, with the number of oscillations increasing as the line 4k + q = is approached. 
At the same time, the SH component keeps a single-humped shape, with an increasing 
amplitude. 

It is interesting to interpret these numerical observations geometrically. Note that 
geometric arguments have been successfully used in Ref. to understand the transi- 
tion from ESs to ordinary gap solitons (however, that work does not explain why this 
transition is absent in the system (|5jl). 

Observe first that in the region k 2 < 1 and 4k + q < 0, where the ESs exist, the 
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O-equilibrium has a two-dimensional center manifold W c , which is filled with periodic 
orbits. In fact, from Eqs. (0) it can be seen that the flow in W c obeys the linear 
equation 

W" = (4k + q)W, 

which gives rise to periodic solutions with arbitrarily large amplitude. This explains 
the (unusual) behavior of the ES solutions, shown in Figure EJ More precisely, for dif- 
ferent branches of the ESs in Eqs. (0), the SH component W grows indefinitely when 
the line 4k + q = is approached. At the same time, the FF component, which corre- 
sponds to the hyperbolic direction of the O-equilibrium, shrinks and develops more and 
more internal oscillations. Thus, in the (Vi, V2, W, W) phase space, the corresponding 
homoclinic orbit follows a periodic orbit in W c , thereby oscillating around W c . Both 
the size of the periodic orbit that is followed and the number of oscillations around 
W c increase as 4k + q = is approached. 

Although these geometric considerations give insight into the behavior of Eqs. ©, 
they do not explain the most remarkable feature, viz., the accumulation of branches of 
the ESs at k = 1, q = —4. In fact, since the O-equilibrium has a four- dimensional center 
manifold for these parameter values, even the local behavior near this equilibrium is 
completely unknown a priori. We therefore address the existence of ES solutions via 
an asymptotic analysis of a normal-form system derived below. 

To this end, we again consider Eqs. (|1J) with c = 0, but this time reduced to the 
invariant subspace U\ — U% = V, U3 = = W. Being interested in the behavior near 
k = 1 and q = —4, we define a = k — l,/3 = 4k + q. Recall that due to the invariance of 
(j3J) with respect to S, the results will not depend on the choice of a particular invariant 
subspace. 

The equations for V and W then read 

-iV - aV - (V - V*) + VW = 0, 

(9) 

w" -pw + \v\ 2 = 0. 

The scaling (5 = ace, £ = (— a)~ 1 / 2 x, V\ = (— a)vi, V 2 = (—a) 3 ^ 2 v 2 , and W = (—a)w 
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transforms this system into 

v[ = (2 + a)v 2 + av 2 w, 



V 2 = V\ + V\W 



(10) 



w" = —ew — v\ + av 2 , 



where the prime now stands for d/dx. We are interested in solutions at a < 1, for 
which \v 2 \ <C \vi\. In this case, Eqs. (ITU|) can be further simplified to 

v[ = 2v 2 , 

v' 2 = Vi + ViW, (11) 
w" = —ew — vf, 

or, setting v = V\, to a system of two second-order equations, 

v" = 2v + 2vw, 



to" = — ew — v 2 . 



(12) 



System (fT2*|) describes stationary solutions in optical models with the x 2 nonlinearity 
of the so-called type I, which means that the system contains only two wave fields, 
rather than three. It has been studied in detail by a number of authors, see the 
™ QQ. The case we a r e — in, with £ > 0, , S g ene rally e,a SSln e d a S tne 
"bright-dark" case, where the latter precisely means that it may have generic solutions 
localized in one direction only. We view Eqs. (fT2"|) as a normal form for a class of 
models supporting embedded solitons through quadratic nonlinearity. 

Remark 3 Note that Eqs. can also be interpreted in the context normal form 
theory, as, for example, introduced in flty - Choosing the coefficient in front of the 
linear term in the equation for v to be another parameter, say fi, instead of the constant 
value 2, Eqs. Mty) describe a normal form for equilibria with fourfold eigenvalue zero 
in the class of reversible systems. It is interesting to note that the general normal form 
contains an additional quadratic term w 2 in the equation forw. The fact that this term 
is missing in hltyl shows the special (degenerate) character of the equations. (Note that 
this also implies a purely linear flow w" = —ew in the centre manifold of the system.) 
Consequent ly, f ormer studies of normal forms near equilibria with fourfold eigenvalue 
zero in \ do not discuss the accumulation of ESs observed in Eqs. lilty). 
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We shall now proceed with an asymptotic analysis to find a new class of solutions to 
the much-studied equations (JT2J); to the best of our knowledge, the numerous previous 
studies of the system did not discuss the solutions that we report on below. Note 
that the solutions have an unfeasibly large w component and were probably missed 
because of that. Furthermore, as stationary solutions to the full type-I x^ model, 
these solitons are likely to be unstable, but they are, nevertheless, of importance for 
the understanding of the normal forms and intrinsic structure of the x models. 

III. ANALYSIS OF EMBEDDED SOLITONS BY MEANS OF THE WKB 
APPROXIMATION 

We are interested in the existence and shape of ES solutions of Eqs. (112)1 in the 
limit e — > +0. After applying the transformation to map solutions of Eqs. (JHJ) into 
those of Eqs. (|12)1 . we find that in this limit the w-component becomes large and wide, 
while the v -component oscillates in an effective potential induced by w(x) within this 
inner core. Note that we focus on detecting even ES solutions of (JT2J), which are the 
fundamental ones. 

A. Inner zone analysis 

We first consider the core of the ES solution, that is, values of x, for which w does 
not decay to 0. Here we we shall look for solutions of Eqs. |T2*|) via the ansatz 

v(x) = v (5x) sin(0(x)), 



Here it is assumed that <5 1, such that v , u> and W2 are slowly varying functions. 

Substitution of the ansatz (|T3~)) into Eqs. ()12j) and equating coefficients in front of 
sin <p and cos <p leads, in the lowest-order approximation, to the following equations: 



w(x) 



w (V5x) + tw 2 (&c) cos(20(x)). 



(13) 




2(wo-l), 



(14) 



d 

dx 





2C 



(15) 
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where C is an arbitrary integration constant, and the factor 2 was introduced for 
convenience. Observe that Eq. (JT5j) requires w > 1. Next, substituting Eqs. (|T3|). 
and (JHJ), (|T5|) into the equation for w in system (JT2J) , and equating coefficients as 
above, yields the following results: 

d2w ° a £ C n 

vl C 

W2 ~ 2(4<P' 2 -e)~ v^l-(8K-l)-£)' 

After the rescaling £ = ^e/dx, Eq. (JTfiJl may be represented as 

d 2 w = dP (w ) , lg , 

i.e., as the equation of motion for a classical particle with mass m = 1 and the coor- 
dinate wo, in a potential 

i 2C 

PW^-W-ij-i^v^i, ^ = v (19) 

In order for Eqs. (jlHj) to describe ES solutions of Eqs. ()12j1 . which vanish at 
infinity, we must equate the conserved energy of the mechanical model (j!8j) to zero, 
(1/2) (u^f + P(w ) = 0. This equation can be solved via the substitution 

™ (£) = l + 77 4 (£), (20) 

which yields an implicit solution, 

P° An 2 

£ = / ; = drj. (21) 

Jo v 72 ^ - 2 V 2 - ^ 

Here, 7/0 is a solution of 

2K - 2rf -rf = 0. (22) 
This implicit solution will be sufficient for the subsequent analysis. 

B. The Bohr-Sommerfeld quantization of the constant C 

Using the above results, we can apply the Bohr-Sommerfeld (BS) quantization 
rule (a part of the semi-classical or WKB approximation in quantum mechanics), 
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which selects a discrete spectrum of values of e, necessary for the existence of the ES 
solution. Recall that the analysis so far has been valid for the case of Wo > 1. If Wq — 1, 
then 0' = 0, because of Eq. (|14jl . and we face a turning-point problem. Moreover, 
wq(x) = 1 leads to a divergence in Eq. (|15jl. which has to be compensated by setting 
sin(0(x)) = 0. Obviously, there are two values of x for which Wq(x) — 1 vanishes. In 
fact, these points describe the transition from the core of the ES to its asymptotic 
tails, where the solution decays to 0. At one of the two points, say x = 0, we simply 
set = 0. Then the condition that sin0 also vanishes at the second point, x = X, 
implies that 



Observe that n counts the number of "inner oscillations" of the f-component of the 
derived solution. We thus obtain 



Equation ()24|) is tantamount to the BS rule in the semi-classical limit of quantum 
mechanics. 

We now have to distinguish two cases. First we assume if < 1. In this case we 




(23) 




(24) 





(25) 



In particular, for small e, making use of the fact that 




(26) 



the BS quantization rule takes an eventual explicit form, 



(27) 



On the other hand, for large K 3> 1, we have 7/q = 2K 



. In a similar way, using the 



numerical value 




(28) 
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the BS rule leads to 



{2Kn ) 1/3 = y/i. (29) 

1 ' Ay/2- 0.373 V 1 ' 



C. The asymptotic solution in the outer zone, and quantization of e 

The above analysis has been dealing with the inner zone (core) of the ES solution, 
where wq took values wq > 1. It remains to match the results to the solution in the 
outer zone (tail), where wq decays to 0. To this end, we note first that the solution 
for the fundamental field, which takes the form of Eq. (jl3|) in the inner zone, can be 
matched to an exponentially decaying solution for the same field in the outer zone, 
given by 

Uouter(z) = Aexp (^-V2\x\j . (30) 

The actual issue is to relate the constant A in this expression to C in Eq. (|15J). This 
can be done using the known connection formula of the WKB approximation, with 
the result 



A = v^72. (31) 

In the same outer zone, we can use the fact that the ES is a localized solution. This 
implies that the corresponding solution for the SH field is 



Sexp (-2V2|ar|) . (32) 



w(x) 



The amplitude B can be computed using the differential equation (jT2jl for w and 
relation (|31jl. We thus find 

B = -A 2 /8 = -C/16. (33) 

Now, we can match the outer and inner solutions for the SH field, demanding 
continuity of this field. According to the above results, the inner-field solution takes 
the value (^) inncr (x = 0) = — wo(x = 0) = —1, see again (JT3J), at the border between 
the two zones. The matching condition thus implies C = 16, and therefore we find 

eK n = 32. (34) 
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It should be noted that this matching procedure is not a rigorous one, which would 
demand to construct an intermediate asymptotic solution in a transient zone, whose 
functional form (rather than just values of the fields at the contact points) must be 
matched to ones in the inner and outer zones. Such a rigorous procedure would be 
extremely complicated in the present system, while the simpler one described above 
yields, eventually, rather accurate results, as shown below. 

We can now use the "quantization rules" ([2~7|) and (J2*9j) to derive approximations 
for values e n at which ES solutions exist in system (JO}. We again start with the case 
K <C 1. The above considerations show that 

32 16^ 



(35) 



and therefore 



However, this implies K n — > oo for n — > oo, in contrast to K <C 1- So this case is 
irrelevant. 

In the opposite case, when K is large, the combination of Eq. ()34|) and the corre- 
sponding approximation (|2T?j) produces the result 

e n = 3.27 -n- 6/5 . (37) 

The applicability of this result also demands n to be large enough. The condition 
K 1, which was used to derive Eq. ()37|). now reduces to 32/e n « 9.8 • n 6//5 ^> 1, 
which is always satisfied. So, Eq. (|3*7j) yields an asymptotic formula for values of e n at 
which Eqs. fjT2"j) support ES solutions. 



D. Numerical verification 



We now investigate the validity of the analytical prediction result (13 7|) through 
numerical solution of Eqs. (|T2"|) . Similar to the case of system (j^J), ES solutions can 
be computed by searching intersections of the unstable manifold of with Fix(R) by 
dint of the shooting method. This way, the first 140 solutions have been constructed. 
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FIG. 4: Enter Figure here. 

The results are summarized in Fig. HJ In the main panel e) we plot n~ 6 / 5 against e, 
where n is identified as the number of oscillations of the w-component in the solution's 
inner zone, following the consideration of the soliton's structure presented above in 
Section ITTT1 For n = 1,2,10,20, the computed solutions are shown in panels a) - d). 
Panel f) specially displays a zoom of the small square from e), which corresponds to 
50 < n < 140. For these values of n we may reasonably assume may Eq. (J37|) to be 
applicable. Again, the axes are chosen to be n -6 / 5 and e with the objective to compare 
the numerical and analytical results. 

A linear regression analysis for the close-up in panel f) yields the following fit to 
the computed values, 

e n = 6 • 10~ 5 + 2.8247 • n" 6/5 , (38) 

which is also shown in Fig. Eff). First of all, it corroborates the main result of the 
above asymptotic analysis, viz., that the scaling index for e n is —6/5. Further, the 
numerical factor in front of n -6 / 5 , predicted by Eq. (J3J), and its empirical counterpart 
in Eq. ()38[) differ by ~ 15%. This moderate discrepancy may be explained by the 
insufficiently accurate matching procedure, as discussed above. Finally, the very small 
constant term in Eq. ()38j) is actually caused by the contribution from the ES solutions 
with smaller n, to which the asymptotic analysis does not apply. In particular, limiting 
the numerical data to 80 < n < 140 reduces the constant term by a factor of 6, and 
shows that all the data fall into a 95%-confidence interval corresponding to a strictly 
proportional relation between e n and n~ 6//5 , without the constant term in Eq. (}3*Bj) . 

IV. CONCLUSIONS 

In this paper, we have derived a system of two second-order equations with one 
free parameter, e, as an asymptotic normal form of models with quadratic nonlinear- 
ity, close to a critical point from which branches of embedded-soliton (ES) solutions 
originate. At the critical point the zero solution possesses a fourfold eigenvalue zero. 
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The derived normal form equations are tantamount to a well-known system, which 
is fundamental for the usual second-harmonic-generation model of type I. Despite the 
fact that the system has been much studied, we were able to find a new infinite series of 
ES solutions in this work. A fundamental characteristic of the series is an asymptotic 
distribution law for discrete values of e at which the ES exist. Making use of the Bohr- 
Sommerfeld quantization rule borrowed from the semi-classical version of quantum 
mechanics and a simplified procedure of matching the solutions in the inner and outer 
zones of the soliton we were able to predict the asymptotic distribution law in a fully 
analytical form. The comparison with direct numerical results shows that the predicted 
scaling law, s n = const • n~ 6//5 , may be absolutely exact (no discrepancy in the value 
of the scaling index, —6/5, was found), and the constant factor was predicted with an 
error ~ 15% (due to an inaccuracy of the simplified matching procedure). 

These results present fundamental (and previously unknown in any form) informa- 
tion about families of ES solutions in multi-wave systems with quadratic nonlinearity. 
Additionally, they shed new light on the above-mentioned system (JT2J) that describes 
type-I second-harmonic generation. 

It is worthwhile to contrast the features of the families of ESs that we have estab- 
lished with other infinite families of ESs found in the literature, e.g. 
all other cases we are aware of the individual members of the family correspond to 
multi- humped solitary waves. In contrast, our construction produces waves that are 
fundamental in one component and become highly oscillatory in the second component. 
Effectively these are single-humped waves. 

There are chiefly two mechanism that produce multi-humped waves, under a so- 
called Birkhoff signature condition on the sign of the nonlinear term. This condition 
is satisfied by the three-wave system investigated here. One theory (see remark 2 
above) shows that there must be an infinite family of branches of multi-humped ESs 
accumulating on each branch of single-humped ES. The other theory ^| shows that 
multi-humped waves must accumulate on the limit corresponding to k — 1 in this 
model. Further numerical results not presented here compute representatives of 
both kinds of multi-humped ESs and follow their branches in the (k, g)-plane. Note 



17 



Accumulation of embedded solitons 



that the latter kind also produces a family of curves that originate from the point 
q = —4, k = 1 in Figure 2 jj, Fig 4.3] similar to the families computed by Yang and 
Akylas for the third-order nonlinear Schrodinger equation. The aim of this paper 
has been to establish why in addition, an infinite family of single-humped solitons 
should emanate from the same critical point in the three-wave model. 

We have already pointed out that we do not expect these newly-established families 
of excited ESs to be stable for the full PDE system. This has been demonstrated 
in numerical studies for the 3- wave system (JJJ in jjj. Nevertheless, it has also been 
shown there that if the parameters are chosen close enough to the critical value, the 
even ground-state solution shown in the top panel of Figure El is virtually stable. This 
means that although the central cores emit radiation tails, those have an extremely 
small amplitude and it takes a considerable distance for the soliton solution to show 
instability. We would expect that the excited states would also develop very slow 
instabilities near to the critical point. Thus, there may be physical application of these 
multitude of states, for example in optical memory, as in most experimental systems 
the residence time is much smaller than that required to produce the instability. 
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List of Captions 

FIG. 1: Results of the linear analysis predicting the possible existence of embedded solitons in 
the case of D = 1. In each part of the parameter plane, the set of eigenvalues found from the 
equations © linearized around the zero solution is displayed. The large dot at {k, q) = (1, —4) 
denotes a codimension-two point at which the linearized system has a fourfold eigenvalue zero. 



G. 2: Branches of even (dotted curves) and odd (solid curves) ESs in system (|5*)). as per Ref. 
4j. No branch crosses the line Ak + q = 0, but all of them approach the point (k, q) = (1, —4) 
tangentially to this line. 



FIG. 3: Examples of even (left panels) and odd (right panels) embedded-soliton solutions 
of Eqs. ((HJ), computed for k = 0, as per Ref. J]. Recall that V\ and V2 are the real and 
imaginary parts of the complex field V. 



FIG. 4: Panels a) through d) show characteristic examples of the embedded solitons found 
from the numerical solution of Eqs. ()12j) . 140 values of e at which the embedded solitons are 
found are collected in the main box e), n being identified as the number of oscillations of the 
^-component in the soliton's inner zone. Panel f) is a close-up of a part of the main diagram 
for 50 < n < 140. It also shows the straight line (|38j) which provides for the best fit to the 
data included in this panel. 
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Figure 1 
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Figure 2 
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Figure 3 
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Figure 4 
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